library(ggplot2)
## Warning: package 'ggplot2' was built under R version 3.5.2
library(corrplot)
## corrplot 0.84 loaded
library(qvalue)
## Warning: package 'qvalue' was built under R version 3.5.2
IBD_ldsc = get(load("/Users/kushaldey/Documents/singlecellLDSC/output/Healthy_gene_score_Feb1_plus_All_IBD.rda"))
dim(IBD_ldsc)
## [1] 51 64 3 4
tau-star
Roadmap Enhancers linked to Genes
tau_table = IBD_ldsc[,,3,1]
ptau_table = IBD_ldsc[,,3,2]
qtau_table = matrix(qvalue(as.vector(ptau_table))$qvalues, nrow = nrow(tau_table), ncol = ncol(tau_table))
ptau_table2 = matrix(qvalue(as.vector(ptau_table))$pvalues, nrow = nrow(tau_table), ncol = ncol(tau_table))
tau_table[which(qtau_table > 0.1)] = 0
#tau_table[tau_table < 0] = 0
corrplot(tau_table, is.corr = F, tl.cex = 0.7, tl.srt = 45)

apply(tau_table, 2, mean)[order(apply(tau_table, 2, mean), decreasing = T)][1:10]
## PASS_IBD_deLange2017
## 0.3114454
## UKB_460K.blood_RBC_DISTRIB_WIDTH
## 0.2301954
## UKB_460K.disease_HYPOTHYROIDISM_SELF_REP
## 0.1990420
## UKB_460K.blood_PLATELET_COUNT
## 0.1922581
## UKB_460K.disease_ALLERGY_ECZEMA_DIAGNOSED
## 0.1760991
## PASS_Rheumatoid_Arthritis
## 0.1670322
## PASS_Alzheimers_Jansen2019
## 0.1575800
## UKB_460K.blood_WHITE_COUNT
## 0.1480575
## UKB_460K.blood_RED_COUNT
## 0.1476459
## UKB_460K.biochemistry_TotalProtein
## 0.1389440
100kb linked to Genes
tau_table = IBD_ldsc[,,3,1]
ptau_table = IBD_ldsc[,,3,2]
qtau_table = matrix(qvalue(as.vector(ptau_table))$qvalues, nrow = nrow(tau_table), ncol = ncol(tau_table))
ptau_table2 = matrix(qvalue(as.vector(ptau_table))$pvalues, nrow = nrow(tau_table), ncol = ncol(tau_table))
tau_table[which(qtau_table > 0.1)] = 0
corrplot(tau_table, is.corr = F, tl.cex = 0.7, tl.srt = 45)

Enrichment
Roadmap Enhancer
E_table = IBD_ldsc[,,3,3]
pE_table = IBD_ldsc[,,3,4]
qE_table = matrix(qvalue(as.vector(pE_table))$qvalues, nrow = nrow(E_table), ncol = ncol(E_table))
pE_table2 = matrix(qvalue(as.vector(pE_table))$pvalues, nrow = nrow(E_table), ncol = ncol(E_table))
E_table[which(qE_table > 0.01 & pE_table2 > 0.001)] = 1
E_table[E_table < 1.1] = 1
Edif_table = E_table - mean(E_table)
Edif_table[Edif_table < 0] = 0
corrplot(Edif_table, is.corr = F, tl.cex = 0.7, tl.srt = 45)

apply(E_table, 2, mean)[order(apply(E_table, 2, mean), decreasing = T)][1:10]
## PASS_IBD_deLange2017
## 12.133796
## UKB_460K.blood_RBC_DISTRIB_WIDTH
## 10.034983
## UKB_460K.biochemistry_TotalProtein
## 10.019738
## UKB_460K.disease_HYPOTHYROIDISM_SELF_REP
## 9.557804
## PASS_LDL
## 9.482519
## UKB_460K.blood_PLATELET_COUNT
## 9.385778
## PASS_Alzheimers_Jansen2019
## 9.143306
## UKB_460K.blood_RED_COUNT
## 8.696251
## PASS_Rheumatoid_Arthritis
## 8.674550
## UKB_460K.blood_WHITE_COUNT
## 8.457306
100kb
E_table = IBD_ldsc[,,1,3]
pE_table = IBD_ldsc[,,1,4]
qE_table = matrix(qvalue(as.vector(pE_table))$qvalues, nrow = nrow(E_table), ncol = ncol(E_table))
pE_table2 = matrix(qvalue(as.vector(pE_table))$pvalues, nrow = nrow(E_table), ncol = ncol(E_table))
E_table[which(qE_table > 0.01 & pE_table2 > 0.001)] = 1
E_table[E_table < 1.1] = 1
Edif_table = E_table - mean(E_table)
Edif_table[Edif_table < 0] = 0
corrplot(Edif_table, is.corr = F, tl.cex = 0.7, tl.srt = 45)

apply(Edif_table, 2, mean)[order(apply(E_table, 2, mean), decreasing = T)][1:10]
## UKB_460K.blood_RBC_DISTRIB_WIDTH
## 0.6817303
## UKB_460K.blood_PLATELET_COUNT
## 0.6014767
## PASS_Alzheimers_Jansen2019
## 0.5986725
## PASS_LDL
## 0.6014308
## UKB_460K.biochemistry_Cholesterol
## 0.5424124
## UKB_460K.biochemistry_TotalProtein
## 0.5272918
## UKB_460K.biochemistry_AlkalinePhosphatase
## 0.5150232
## UKB_460K.biochemistry_TotalBilirubin
## 0.4795640
## PASS_IBD_deLange2017
## 0.4669224
## UKB_460K.biochemistry_AspartateAminotransferase
## 0.4445825